{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "np.random.seed(666)\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.empty((100,2))\n",
    "X[:,0] = np.random.uniform(0., 100., size=100)\n",
    "X[:,1] = 0.75 * X[:,0] + 3. + np.random.normal(0., 10., size=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def demean(X):\n",
    "    return X - np.mean(X, axis=0)\n",
    "\n",
    "X = demean(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fa6a04c2550>]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGbCAYAAAAY8u5bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df4yd1X3n8c+XsSGTXzskNgFf49iRiLWhTjCZsuxaq3oJ1ISkYDlplt1uy6bRWluRVVilDnaQtkRq5WnZhiRqS9dKIiVStgQ11HELXRfieKVEC8mYgVCHurFCCx6zAboMbeqpsc13/7jPxdczz5374znPj/M875eEPPfHzD1zH8z9cM73fI+5uwAAABDOeWUPAAAAoG4IWAAAAIERsAAAAAIjYAEAAARGwAIAAAhsWdkD6LZixQpfu3Zt2cMAAADo69ChQy+6+8q0xyoVsNauXavp6emyhwEAANCXmf1tr8dYIgQAAAiMgAUAABAYAQsAACAwAhYAAEBgBCwAAIDACFgAAACBEbAAAAACI2ABAAAERsACAAAIjIAFAAAQGAELAAAgMAIWAABAYAQsAACAwJaVPQAAAHrZOzOru/Yf0fG5ea2aGNeOLeu1dWOr7GEBfRGwAACVtHdmVrvuf1Lzp85Ikmbn5rXr/icliZCFymOJEEA93XZb+x9E6679R14LVx3zp87orv1HShoRMDhmsADU0+OPlz0CZHR8bn6o+4EqYQYLAFBJqybGh7ofqBICFgCgknZsWa/x5WPn3De+fEw7tqwvaUTA4FgiBABUUqeQnV2EiBEBCwBQWVs3tghUiBJLhAAAAIERsAAAAAIjYAEAAARGwAIAAAiMgAUAABAYuwgBAGgoDtPODwELAIAG4jDtfLFECABAA3GYdr4IWAAANBCHaeeLgAUAQANxmHa+CFgAADRQbIdp752Z1aapA1q38wFtmjqgvTOzZQ9pSRS5AwDQQDEdph1jQT4BCwCAhiryMO0sLSGWKsgnYAEAgEbKOgMVY0E+NVgAAEQgthqkbllbQsRYkE/AAgCg4jozQLNz83KdnQGKJWRlnYGKrSBfChiwzGzMzGbM7M+S2+vM7FEz+5GZfd3Mzg/1WgAANEnsTUGzzkBt3djS7m0b1JoYl0lqTYxr97YNla2/ksLWYH1C0lOS3pzc/m1Jd7v7vWb2h5I+JumegK8HAEAjxFiD1G3HlvXn1GBJw89ADVqQX5XzFYPMYJnZakkfkPTF5LZJukbSHydP+YqkrSFeCwCApomxBqlbUTNQVVpKDTWD9TlJn5L0puT2WyXNufvp5PYxSanvopltl7RdktasWRNoOAAA1EeIGaCyFdESokrtHDLPYJnZByU97+6Huu9Oeaqnfb+773H3SXefXLlyZdbhAABQOzHWIJWhSkupIWawNkm60cxukPQ6tWuwPidpwsyWJbNYqyUdD/BaAAA0UpFNQWO1amJcsylhqoyl1MwzWO6+y91Xu/taSTdLOuDuvyTp25I+nDztFknfzPpaAAAAvVSpnUOendxvl3Svmf2mpBlJX8rxtQAAQEZV2YE3qiqdr2juqaVRpZicnPTp6emyhwGgDjZvbv958GCZowCisfA4G0laPmZ6w/nL9PL8qSgDV97M7JC7T6Y9xlmEAAAgdQfeqTOuuflTkoY/P7DpOCoHAAAMtNMupu7xZSNgAQCAgXfaxdI9vmwELAAAkLoDL00s3ePLRsACAACLmple+PrlWn7euX3DY+seXyaK3AEAgKTFzUxjb9tQJgIWAABIRff40bFECAAAEBgzWAAAoFLqsDRJwAIAAJWxsKN8rA1OWSIEAACVkdZRPsYGpwQsAABQGb0amcbW4JSABQAAKqNXI9PYGpwSsAAAaKi9M7PaNHVA63Y+oE1TB7R3ZrbsIaV2lI+xwSlF7gAANFDVism7dw5OvH65Llh2nl6eP8UuQgAAUI5R2hosVUxedJhZGPZeOnFK48vHdPe/vSK6YNXBEiEAABHrhJPZuXm5zs5E9Vvuq0ox+d6ZWX3yvidqsXOwGwELAICIjdrWoArF5J1weMY99fHYdg52I2ABABCxUWeiqlBMnhYOu8W2c7AbAQsAgIiNOhO1dWNLu7dtUGtiXCapNTGu3ds2FFrztFQIjHHnYDeK3AEAiNiOLevPKRCXBg8nWze2Si0iXzUxrtmUkDVmVnjYC40ZLAAAIlaFmahR9Vqm/N2PvCeK8S+FGSwAACJX9kzUqDpjHrbFRAwIWAAAoDSxhsN+WCIEAAAIjIAFAAAQGAELAAAgMGqwACCjUc6BA1BvBCwAyGDhIbWdc+AkEbKABmOJEAAyGPUcOAD1xgwWAGQw6jlwqJYmLfM26XctEwELADLoddRHzIfUNk1Tlnn3zszqM396WC+dOPXafXX9XauAJUIAyKDXUR8xH1LbNL2WeT953xPaOzNb0qjC6oTI7nDVUfSS9t6ZWW2aOqB1Ox/QpqkDtXmPF2IGCwAyqPNRH1WUx/JWr+XcM+6lzu6E/F3TQmS3opa0mzJbKBGwACCzuh71UTV5fTj3WuaVzp3dKTJEh/5d+wWoopa0l9oUUre/QywRAgCikNeOzbRl3m6dcDM7Ny/vup3n0lbo33WpAFXkknaTNoUQsAAAUcjrw3nrxpZ2b9ugMbPUx8fMCm/FEfp37RUiJ8aXa/e2DYXNHvUKenXcFMISIQDgNVXewp/njs3O79i9LCe1Z3d61S7lOesS+netSq3gji3rU9/jOm4KIWABACQVW4DcL8ilPZ73h3OvEHLX/iOFt+LI43etQq1gVYJeEQhYAABJxRUg9wtyvR7fvW2DPvTelv7o0Wd1xl1jZvrQe8OGhl4hpOhZlzoHkSoEvSIQsAAAkrLV/QyztNgvyPV6/M59h3Xy9Ks64y6p3UbhG4dmNfn2t+T6gV1W2GlKEKkrAhYAQNLodT/DLi32C3K9Hp+b790kk7CDqmEXIQBA0uhd6YdtKdBvJ9mwtU113OKP+BGwAKDmBj2apNOuoDUxLpPUmhgfaAv/sEuL/YJcr8cvfP3y1J9Xxy3+GFxVj97JvERoZpdK+qqkiyW9KmmPu3/ezN4i6euS1kr6G0kfcfeXsr4eAGBwwy7fjbIUNuzSYr+apl6PS8UXm6Paqnz0TogarNOSPunuj5nZmyQdMrOHJP1HSd9y9ykz2ylpp6TbA7weAGBARewMHKWlQL8gt9TjddxZh9FU+eidzAHL3Z+T9Fzy9T+Y2VOSWpJukrQ5edpXJB0UAQsAClXE0SRF7rKj2Bzdqnz0TtBdhGa2VtJGSY9KelsSvuTuz5nZRT2+Z7uk7ZK0Zs2akMMBgMbLs/t5t7oEnyp3ssdiRf37PYpgRe5m9kZJ35B0m7v//aDf5+573H3S3SdXrlwZajgAAI2+MzBmoxY9d+p5ijzUGdlU+d/vIAHLzJarHa6+5u73J3f/xMwuSR6/RNLzIV4LADCc1y0/+5/6og/3LVqWkDRsuwmUb9Sdr0UIsYvQJH1J0lPu/tmuh/ZJukXSVPLnN7O+FgBgcAt3WEnSydOvljii/N257/DIRc951vOw9Jifqi5Ph5jB2iTplyVdY2aPJ//coHawus7MfiTpuuQ2AKAgTZuR2Tszm9rtXRosJPVrgJplXCw9Nk+IXYTfkWQ9Hn5f1p8PABhNlXdY5WGp4DhISBql3cSg4+p1tiKzWvXFWYQAUFNV3mGVh6WC4yAhKa92E0udrdiZcatSg0yEQcACgJrKa0amqnoFygtfv3zg0JJHPU+vcS1UlQaZCIOzCAGgpqq8wyoPvbbs/8YvXF7SiNrSxtVLXZdvm4gZLACosarusMpDkR3ls47rxCun9dKJxQX5dV2+bSICFgCgNqoaKBeOK62FRp2Xb5uIgAUAQJcielZVdbYN4RCwAABILJxZynN3X1Vn2xAGRe4AACSa1pwV+SFgAQCQaFpzVuSHgAUAQCKv43LQPAQsAEAp9s7MatPUAa3b+YA2TR2oxNl8vXppsbsPw6LIHQBQuCKLyYfB7j6EQsACkEkRW9pRP0sVk5f97w+7+xACAQvAyKo6C4Hqo5gcdUcNFoCRsaUdowpdTF7Fei40GwELwMjqPAvBB3a+QhaTd2ZSZ+fm5To7k8o1Q5kIWABGVtct7Xxg52/rxpZ2b9ug1sS4TFJrYly7t20YaWmZmVRUETVYAEa2Y8v6Wh5YW+UC7DoJVUxe55lUxIuABWBkdd3Szgd2dQyyS3XVxLhmU65N7DOpiBsBC0AmddzSzgd2NQy6S7WuM6mIGzVYALAA3byrYdDaqpD1XEAozGABwAJ1XfqMzTBLtXWcSUXcCFgAkIIP7PL1Wqo9z0zrdj5A8EWlsUQIAKiktKVaSTrjTvsMVB4BCwBQSQtrq8bMFj2HfleoKpYIAQCV1b1Uu27nA6nPoX0GqogZLABAFOp6cgDqiYAFAIgC7TMQE5YIAQBRoH0GYkLAAgBEg/YZiAUBC0Bmg5wXBwBNQsACkMmg58UBQJNQ5A4gk0HPiwOAJiFgAchkmPPiAKApCFgAMqE3EQAsRsACkAm9iQBgMYrcAWRCbyIAWIyABSCzqvcmoo0EgKIRsADUGm0kqovgizojYAGotaXaSPBhnt2oIYngi7ojYAGotbLbSNR5liZLSCL4ou7YRQig1spsI9EJILNz83KdDSB7Z2Zzf+0iZGkyW3bwBfJGwAJQa2W2kah7l/ssIYn+aag7AhYab+/MrDZNHdC6nQ9o09SB2swuoG3rxpZ2b9ug1sS4TFJrYly7t20oZBmq7rM0WUIS/dNQd9RgodEotK2+EDVMZbWRWDUxrtmUMFWXWZodW9af8/dHGjwk0T8NdZd7wDKz6yV9XtKYpC+6+1TerwkMikLbaos9AGcJIDHIGpKq3j8NyCLXgGVmY5J+X9J1ko5J+r6Z7XP3H+b5usCgYlnCqfNOtKXEHoCbMEtDSALSmbvn98PN/qWkO919S3J7lyS5++60509OTvr09HRu45Ek/eQ26Z8ez/c1EI3Hnp3TK6fPLLr//GVjuvLSiRJGtNiLPz2pH7/4j3q16+/qeWZ6x4o3aMUbLyhxZPl75Om/6/nY1eveuvQ3P5H8PX/PFQFHBKASXneF9LbPlT0Kmdkhd59MeyzvIveWpGe7bh9L7nuNmW03s2kzm37hhRdyHg5wrjUXjus8s3PuO89May6sTo3MMy/NnxOuJOlVdz3zUrVm2dK8+NOTeuzZOT3y9N/psWfn9OJPTw71/ecvGxvqfgCoirxrsCzlvnM+Kdx9j6Q9UnsGK+fxVCLxojpWSPpOyvLbVZdXZ8njQ/c8oLS/GCbp6akPFD2cgS2sn5La9UfD7OB75v/1/hlXvr3Pz7hlc/vPgweHHDkAZJd3wDom6dKu26slHc/5NYGhVL2GJNadaCHqp5pQwwSgnvIOWN+XdJmZrZM0K+lmSf8+59cEaiXWnWihNhBUPQCH1tQNDUDd5Bqw3P20mX1c0n612zR82d0P5/maQN1UdRanXxAYZOaNMHGu2NtSADgr9z5Y7v6gpAfzfh2gn5g/zKs2izNIEOg380aYWCz2thQAzuKoHDRCkYfuNuHonUHO2Ot3RE3dz+kbRSx92QD0x1E5aISiZgaaMiszaBBYauYttjBRxAxorBsaACzGDBYaoagP86bMymQ55DfkzyhKUTOgHIAM1AcBC41Q1Id5bLMyowoRBGIKE0UF537LqgDiwRIhopB1eaaoVgd1WOIZ5L0OsbOx7N2Rw/w7VWRwrtqGBgCjIWChEpb6sAtR11TUh3msPas6hnmvQwSBssLEsP9O1SE4AygWAQul6/dhF6pAvYgP87JnZbJqSpuAYX/P2IMzgOIRsFC6fh92sdU1xbzEE9t7Paphf8/YgzOA4hGwULp+H3YszxSnKe/1KL9nzMEZQPHYRYjS9dvhF9Nus9g15b1uyu8JoDzMYCGYUXf69atvYXkmm2GuS1Pe66b8ngDKY+5e9hheMzk56dPT02UPAyNYWKgutUPSoD18Yj4nsMqyXpeobd7c/vPgwTJHAaDGzOyQu0+mPcYMFoLIuvuM+pZ8NGVXIABUDQELQTRl91lsyrwuzEoCaDKK3BFETOfKNUlZ16Wos/sAoKoIWAiCXVnVVNZ1acqh1wDQC0uECKLfriyWi8pR1m45lowBNB0BC8H0KlQPcZYgRlfGBoKmNCwFgF5YIkTuWC5qHpaMATQdM1jIHctFzUMjTwBNR8BC7lguaqYYeptRGwggLywRIncsF6GKaCUBIE8ELORu68aWdm/boNbEuExSa2K8GUe1oNKoDQSQJ5YIUYgYlovQLNQGAsgTM1gAGonTBwDkiYAFoJGoDQSQJ5YIATQSrSQA5ImABaCxqA0EkBcCFoDS0Y8KQN0QsACUirMqAdQRRe4ASkU/KgB1RMACUCr6UQGoI5YIAWSStX6KsyoB1BEzWABGFuI8P/pRAagjAhaAkYWon+KsSgB1xBKh2CIOjCpU/RT9qADUTeNnsEIscQBNxXl+AJCu8QGLLeLA6KifAoB0jV8iZIs4MDrO8wOAdI0PWGwRB7KhfgoAFmv8EiFLHAAAILTGz2DVaYmD3ZAAAFRD4wOWVI8lDg7MBQCgOhq/RFgX7IYEAKA6mMGqCXZDNhdLwwBQPZlmsMzsLjP7KzP7gZn9iZlNdD22y8yOmtkRM9uSfahYCg0fm4lGuQBQTVmXCB+S9DPu/m5Jfy1plySZ2bsk3SzpcknXS/oDMxvr+VOQGbsh47d3Zlabpg5o3c4HtGnqwEAhiaXhc3W/h489M6cXf3qy7CEBaKhMS4Tu/hddNx+R9OHk65sk3evuJyU9bWZHJV0l6f9keT30VqfdkFnFuGQ26iYFlobPWvgevnL6jH78wj/qOzOzlb/+AOonZA3Wr0r6evJ1S+3A1XEsuW8RM9suabskrVmzJuBwmqcOuyGzinU35VIzUUuNm0a5Z6W9h6+6930PASAPfZcIzexhM/vLlH9u6nrOHZJOS/pa566UH+VpP9/d97j7pLtPrly5cpTfAQ2z1FJarEtmo85EsTR8FrN5AKqk7wyWu1+71ONmdoukD0p6n7t3QtQxSZd2PW21pOOjDhLo6DdDFeuH7KgzUSwNn8VsHoAqybREaGbXS7pd0s+5+4muh/ZJ+p9m9llJqyRdJul7WV6riWKsJcpbv6W0sj5ks16rHVvWnxMcpcFnolgabkt7D88za+RsHoDyZd1F+HuS3iTpITN73Mz+UJLc/bCk+yT9UNL/knSru5/p/WOwENvv0/WboSpjySzEtdq6saXd2zaoNTEuk9SaGNfubRsITkNY+B6ev2xM71j5Bt5DAKWws6t65ZucnPTp6emyh1EJm6YOpM7EtCbG9d2d15QwomoY5H0peuaPa1VRmze3/zx4sMxRAKgxMzvk7pNpj9HJvaJirSXK2yBLaUUvmXGtAAALcRZhRdGZPV0Vl9K4VgCAhZjBCiCPJaksRc91V7Wibq4VAGAhAlZGeTW2ZPt9MUKEY64VAGAhitwzosA5XgvDsdSeeSp7yRGBUOQOIGdLFblTg5URBc7xirXrOwCg+ghYGVHgHC/CMQAgLwSsjJp2FtxS5wDGhnAMAMgLASujKrYNyEvduss3LRwDAIpDkTsG1qugf8xMr7pHuXuO8x5rjCJ3ADmjkzuC6FWbdCYJ6aFaVBSpaj21kB/CNIAisUSIgQ1Sm8QuPFRR3Za3AVQfAQsDS6tZSjPqLrw6FdCjWmjJAaBoLBHmqG5LEgs7lp9n9tryYLdRduHl1REfkGjJAaB4zGDlpK5LEls3tvTdndfo6akP6Hc/8p5gu/CYYUCeaMkBoGgErJz0Cgx37jtc0ojCC9mighkG5ImWHACKxhJhTnoFg7n5U9o7M1ubZa9Qu/BWTYyntoBghgEhcCA3gKIRsHLSKzBI7f/I8x/2c+3Ysj714GVmGBAKLTkAFIklwpwsFQxY9lqsSR3xAQD1xwxWTrZubOkzf3pYL504teixpZa96rbzcBjMMAAA6oIZrBz9xi9cPlRhbV13HjYBPbwAAN0IWDkadtmLVgVxIhgDABZiiTBnwyx70aogTksFY5Y8AaCZmMGqEJohxolgDABYiIBVIU1uhhhzDRPBGACwEAGrQpraqiD2GqYmB2MAQDpqsCqmia0KYq9hoks4AGAhAhZKV4capiYGYwBAbywRonTUMAEA6oaApbgLrOuAGiYAQN00fomwU2DdqQHqFFhLin7JJ5Zjd6hhAgDUTeMDVuwF1r3EFhyrWMMUS0AFAFRP45cI61BgnYZjd7KJvXUEAKBcjQ9YdS2wrmtwLAoBFQCQReMD1qAF1rEVwvcKiC5FMf6yEVABAFk0PmAN0j09xuWitODYEcP4y1bXmU0AQDEaX+QuLS6w7sxWdYqbT7xyOrpC+O6debMpsy5VH3/ZdmxZf84mAYnWEQCAwTV+BmuhtNmql06cSn1u1ZeLtm5s6bs7r5H1eLzq4y9TU8+FBACEwQzWAmnFzb3Esly0amI8dRYrlvGXpYqtIwAAcWAGa4FBZ3ViWi6iUzoAAMUiYC3Qa1ZnYnx5tMtFLHcBAFAslggX6FXcfOeNl1c2kAzScbyKy110SgcA1BUBa4HYzsWL7UicjljHDQDAIMzdyx7DayYnJ316errsYURl09SB1AL2ifHlesMFyyobEnuNuzUxru/uvKaEEaF2Nm9u/3nwYJmjAFBjZnbI3SfTHgtSg2Vmv25mbmYrkttmZl8ws6Nm9gMzuzLE62CxXkX5c/OnKt0YlU7pAIA6yxywzOxSSddJeqbr7vdLuiz5Z7uke7K+DtIN2mqhaufo0SkdAFBnIWaw7pb0KbWPueu4SdJXve0RSRNmdkmA18ICSx2Js1CVZodoHQEAqLNMRe5mdqOkWXd/wuycfuEtSc923T6W3Pdcys/YrvYsl9asWZNlOI2UVpR/4pXTqd3nqzQ7FNtmAgAAhtE3YJnZw5IuTnnoDkmflvTzad+Wcl9qNb2775G0R2oXufcbDxZLO0sxhnP0qtg6AgCAEPoGLHe/Nu1+M9sgaZ2kzuzVakmPmdlVas9YXdr19NWSjmceLQbC7BAAAOUaeYnQ3Z+UdFHntpn9jaRJd3/RzPZJ+riZ3SvpX0h62d0XLQ8iP7HMDtFsFABQR3k1Gn1Q0g2Sjko6IemjOb0OIkazUQBAXQULWO6+tutrl3RrqJ+Nerpr/5Fz6sSks+0kCFgAgJhx2DNKQ7NRAEBdEbBQGpqNAgDqioCF0tBsFABQV3kVuQN90U4CAFBXBCyUKpZ2EgAADIMlQgAAgMAIWAAAAIERsAAAAAIjYAEAAARGwAIAAAiMgAUAABAYbRqWsHdmlh5NAABgaASsHvbOzGrX/U++dhjx7Ny8dt3/pCQRsgAAwJJYIuzhrv1HXgtXHfOnzuiu/UdKGhEAAIgFAauH43PzQ90PAADQQcDqYdXE+FD3AwAAdBCwetixZb3Gl4+dc9/48jHt2LK+pBEBAIBYUOTeQ6eQnV2EAABgWASsJWzd2CJQAQCAobFECAAAEBgBCwAAIDCWCJdAJ3cAADAKAlYPdHIHAACjYomwBzq5AwCAURGweqCTOwAAGBUBqwc6uQMAgFERsHqgkzsAABgVRe490MkdAACMioC1BDq5AwCAUbBECAAAEBgBCwAAIDACFgAAQGAELAAAgMAIWAAAAIERsAAAAAIjYAEAAARGwAIAAAiMgAUAABAYAQsAACAwAhYAAEBgBCwAAIDACFgAAACBEbAAAAACI2ABAAAERsACAAAILHPAMrP/YmZHzOywmf1O1/27zOxo8tiWrK8DAAAQi2VZvtnM/o2kmyS9291PmtlFyf3vknSzpMslrZL0sJm9093PZB0wAABA1WWdwfo1SVPuflKS3P355P6bJN3r7ifd/WlJRyVdlfG1AAAAopA1YL1T0r82s0fN7H+b2c8m97ckPdv1vGPJfYuY2XYzmzaz6RdeeCHjcAAAAMrXd4nQzB6WdHHKQ3ck33+hpKsl/ayk+8zsHZIs5fme9vPdfY+kPZI0OTmZ+hwAAICY9A1Y7n5tr8fM7Nck3e/uLul7ZvaqpBVqz1hd2vXU1ZKOZxwrAABAFLIuEe6VdI0kmdk7JZ0v6UVJ+yTdbGYXmNk6SZdJ+l7G1wIAAIhCpl2Ekr4s6ctm9peSXpF0SzKbddjM7pP0Q0mnJd3KDkIAANAUmQKWu78i6T/0eOy3JP1Wlp8PAAAQIzq5AwAABEbAAgAACIyABQAAEBgBCwAAIDACFgAAQGAELAAAgMAIWAAAAIERsAAAAAIjYAEAAARGwAIAAAiMgAUAABAYAQsAACAwAhYAAEBgBCwAAIDACFgAAACBEbAAAAACI2ABAAAERsACAAAIjIAFAAAQGAELAAAgMAIWAABAYAQsAACAwAhYAAAAgRGwAAAAAiNgAQAABEbAAgAACIyABQAAEBgBCwAAIDACFgAAQGAELAAAgMAIWAAAAIERsAAAAAIjYAEAAARGwAIAAAiMgAUAABAYAQsAACAwAhYAAEBgBCwAAIDACFgAAACBEbAAAAACI2ABAAAERsACAAAIjIAFAAAQGAELAAAgsEwBy8yuMLNHzOxxM5s2s6uS+83MvmBmR83sB2Z2ZZjhAgAAVF/WGazfkfQZd79C0n9LbkvS+yVdlvyzXdI9GV8HAAAgGlkDlkt6c/L1P5N0PPn6Jklf9bZHJE2Y2SUZXwsAACAKyzJ+/22S9pvZf1c7rP2r5P6WpGe7nncsue+5jK8HAABQeX0Dlpk9LOnilIfukPQ+Sf/V3b9hZh+R9CVJ10qylOd7j5+/Xe1lRK1Zs2bAYQNAH1dcUfYIADSYuafmnsG+2exlSRPu7mZmkl529zeb2f+QdNDd/yh53hFJm919yRmsyclJn56eHnk8AAAARTGzQ+4+mfZY1hqs45J+Lvn6Gkk/Sr7eJ+lXkt2EV6sdvFgeBAAAjZC1Bus/Sfq8mS2T9E9KlvokPSjpBklHJZ2Q9NGMrwMAABCNTAHL3b8j6b0p97ukW7P8bAAAgFjRyR0AACAwAhYAAEBgBCwAAIDACLX2EH8AAARhSURBVFgAAACBEbAAAAACI2ABAAAERsACAAAIjIAFAAAQGAELAAAgMAIWAABAYAQsAACAwAhYAAAAgVn7XOZqMLMXJP1t2eOIyApJL5Y9CJyDa1ItXI/q4ZpUD9dkdG9395VpD1QqYGE4Zjbt7pNljwNncU2qhetRPVyT6uGa5IMlQgAAgMAIWAAAAIERsOK2p+wBYBGuSbVwPaqHa1I9XJMcUIMFAAAQGDNYAAAAgRGwAAAAAiNgRczMft3M3MxWJLfNzL5gZkfN7AdmdmXZY2wCM7vLzP4qec//xMwmuh7blVyPI2a2pcxxNo2ZXZ+870fNbGfZ42kiM7vUzL5tZk+Z2WEz+0Ry/1vM7CEz+1Hy54Vlj7VJzGzMzGbM7M+S2+vM7NHkenzdzM4ve4x1QMCKlJldKuk6Sc903f1+SZcl/2yXdE8JQ2uihyT9jLu/W9JfS9olSWb2Lkk3S7pc0vWS/sDMxkobZYMk7/Pvq/134l2S/l1yPVCs05I+6e7/XNLVkm5NrsNOSd9y98skfSu5jeJ8QtJTXbd/W9LdyfV4SdLHShlVzRCw4nW3pE9J6t6lcJOkr3rbI5ImzOySUkbXIO7+F+5+Orn5iKTVydc3SbrX3U+6+9OSjkq6qowxNtBVko66+4/d/RVJ96p9PVAgd3/O3R9Lvv4HtT/UW2pfi68kT/uKpK3ljLB5zGy1pA9I+mJy2yRdI+mPk6dwPQIhYEXIzG6UNOvuTyx4qCXp2a7bx5L7UJxflfTnyddcj/Lw3leMma2VtFHSo5Le5u7PSe0QJumi8kbWOJ9T+3/OX01uv1XSXNf/JPJ3JZBlZQ8A6czsYUkXpzx0h6RPS/r5tG9LuY8+HAEsdT3c/ZvJc+5Qe0nka51vS3k+16MYvPcVYmZvlPQNSbe5+9+3J01QNDP7oKTn3f2QmW3u3J3yVP6uBEDAqih3vzbtfjPbIGmdpCeS/0itlvSYmV2l9v95XNr19NWSjuc81EbodT06zOwWSR+U9D4/21yO61Ee3vuKMLPlaoerr7n7/cndPzGzS9z9uaSM4fnyRtgomyTdaGY3SHqdpDerPaM1YWbLklks/q4EwhJhZNz9SXe/yN3XuvtatT9IrnT3/ytpn6RfSXYTXi3p5c40PPJjZtdLul3Sje5+ouuhfZJuNrMLzGyd2psPvlfGGBvo+5IuS3ZHna/2ZoN9JY+pcZL6ni9JesrdP9v10D5JtyRf3yLpm0WPrYncfZe7r04+O26WdMDdf0nStyV9OHka1yMQZrDq5UFJN6hdTH1C0kfLHU5j/J6kCyQ9lMwqPuLu/9ndD5vZfZJ+qPbS4a3ufqbEcTaGu582s49L2i9pTNKX3f1wycNqok2SflnSk2b2eHLfpyVNSbrPzD6m9k7oXyxpfGi7XdK9ZvabkmbUDsXIiKNyAAAAAmOJEAAAIDACFgAAQGAELAAAgMAIWAAAAIERsAAAAAIjYAEAAARGwAIAAAjs/wODiDHGW3ePQwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,7))\n",
    "plt.scatter(X[:,0],X[:,1])\n",
    "plt.plot([0 for i in np.empty((130,1))],range(-80,50), color=\"r\")\n",
    "plt.plot(range(-50,50),[0 for i in np.empty((100,1))], color=\"gold\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(w, X):\n",
    "    return np.sum((X.dot(w)**2)) / len(X)\n",
    "\n",
    "def df(w, X):\n",
    "    return X.T.dot(X.dot(w)) * 2. / len(X)\n",
    "\n",
    "def direction(w):\n",
    "    return w / np.linalg.norm(w)\n",
    "\n",
    "# 求出第一主成分\n",
    "def first_component(X, initial_w, eta, n_iters=1e4, epsilon=1e-8):\n",
    "    \n",
    "    w = direction(initial_w)\n",
    "    cur_iter = 0\n",
    "    \n",
    "    while cur_iter < n_iters:\n",
    "        gradient = df(w, X)\n",
    "        last_w = w\n",
    "        w = w + eta * gradient\n",
    "        w = direction(w)\n",
    "        if(abs(f(w, X) - f(last_w, X)) < epsilon):\n",
    "            break\n",
    "        \n",
    "        cur_iter += 1\n",
    "        \n",
    "    return w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.77660914, 0.62998273])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "initial_w = np.random.random(X.shape[1])\n",
    "eta = 0.01\n",
    "w = first_component(X, initial_w, eta)\n",
    "w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 去掉第一主成分"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "X2 = np.empty(X.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for i in range(len(X)):\n",
    "#     X2[i] = X[i] - X[i].dot(w) * w\n",
    "# 上式可以进一步向量化，等价于\n",
    "X2 = X - X.dot(w).reshape(-1,1) * w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD7CAYAAABt0P8jAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAYWklEQVR4nO3df2wc5Z3H8c83xnBL2pODMDTZxpcUpT6Rc7GvFhBFrYDSmhb1WCKlIYIT1VWklcrpQpF1TosK5dLGqkWJdNdrCyoCiRSStsFwDcL8SFvUiLR1zg5OylkkQEI2EXEv+FrBCozzvT+86y7Orr3rndlf835JK+/Ozu48w7KfnXznmecxdxcAIDoWVLoBAIDyIvgBIGIIfgCIGIIfACKG4AeAiCH4ASBiSg5+M1tqZr80s5fM7KCZ/Ut6+Xlm9oyZvZz+u6j05gIASmWl9uM3s8WSFrv7f5vZByXtk5SQ9EVJp9y918x6JC1y938ttcEAgNKUHPxnvKHZ45L+I327wt1PpH8cfuXurbO99vzzz/dly5YF2h4AqHf79u37o7s3F7r+WUFu3MyWSeqQ9FtJF7r7CUlKh/8FeV6zQdIGSWppadHg4GCQTQKAumdmR4pZP7CTu2b2AUk/l7TR3f9U6Ovc/T5373T3zubmgn+wAADzFEjwm1mjpkJ/m7vvTC9+I13iyZwHOBnEtgAApQmiV49J+rGkl9z9e1lPPSHp5vT9myU9Xuq2AAClC6LGv1rSP0oaMbPh9LKvS+qVtMPMviTpqKS1AWwLAFCikoPf3X8jyfI8/alS3x8AEKxAe/VUSv9QUn0Dozo+ntKSppi6u1qV6IhXulkAUJVqPvj7h5LatHNEqYlJSVJyPKVNO0ckifAHgBxqfqyevoHR6dDPSE1Mqm9gtEItAoDqVvPBf3w8VdRyAIi6mg/+JU2xopYDQNTVfPB3d7Uq1tjwvmWxxgZ1d806LBAARFbNn9zNnMClVw8AFKbmg1+aCn+CHgAKU/OlHgBAcQh+AIgYgh8AIobgB4CIIfgBIGIIfgCIGIIfACKG4AeAiCH4ASBiCH4AiBiCHwAiJpDgN7MHzOykmR3IWnaXmSXNbDh9+1wQ2wIAlCaoI/4HJV2TY/m97t6evj0Z0LYAACUIJPjd/XlJp4J4LwBAuMKu8d9qZi+mS0GLcq1gZhvMbNDMBsfGxkJuDgAgzOD/gaSLJLVLOiHpnlwruft97t7p7p3Nzc0hNgcAIIUY/O7+hrtPuvtpSfdLujSsbQEAChfaDFxmttjdT6QfXi/pwGzr15v+oSTTQQKoSoEEv5k9IukKSeeb2TFJd0q6wszaJbmk1yR9OYht1YL+oaQ27RxRamJSkpQcT2nTzhFJIvwBVFwgwe/u63Ms/nEQ712L+gZGp0M/IzUxqb6BUYIfQMVx5W4Ijo+niloOAOUUWo0/ypY0xZTMEfJLmmLU/gFUHEf8IejualWsseF9y2KNDbryb5u1aeeIkuMpuf5S++8fSlamoQAiieAPQaIjri1r2hRviskkxZti2rKmTb/8n7G8tX8AKBdKPSFJdMTPKOHctn0457rU/gGUE0f8ZbSkKVbUcgAIA8FfRvlq/91drRVqEYAootRTRpnSD716AFQSwV9muWr/AFBOBH+Vo98/gKCZu1e6DdM6Ozt9cHCw0s2oGjPH/JEk09TgR3F+BACkmdk+d+8sdH1O7laxXGP+ZH6mufgLwHwR/FVsrv79qYlJbdw+rNW9u/kBAFAwgr+KFdq/Pzme0sbtw+q4+2l+AADMieCvYrn6/c/mzbcnKP8AmBPBX8Wyx/yRpk7szoXyD4C50J2zymX3+8907cw15PNMzPoFIB+O+GtIoiOuPT1Xaeu69oJKQIz8CSAXgr8GZUpATbHGOddNjqe0vGcXpR8A0wIJfjN7wMxOmtmBrGXnmdkzZvZy+u+iILaFKYmOuIbv/Iy2rmufPgeQD5O+AMgW1BH/g5KumbGsR9Jz7r5C0nPpxwhYMeUfSj8ApICC392fl3RqxuLrJD2Uvv+QpEQQ20JuM2f9yodJXwAENlaPmS2T9At3/7v043F3b8p6/k13P6PcY2YbJG2QpJaWlo8fOXIkkPZE3ere3Tl7/zTFGrXwnLMY9A2oIzU3Vo+73+fune7e2dzcXOnm1I1cF381LjC99e57TPYORFyYwf+GmS2WpPTfkyFuCzPkmvD9A391liYm3/8vPOr+QPSEeQHXE5JultSb/vt4iNtCDjMnfVnesyvnetT9gWgJqjvnI5JekNRqZsfM7EuaCvxPm9nLkj6dfowKYrJ3AFJAR/zuvj7PU58K4v0RjO6u1jMmdmGydyB6GKsnQpjsHYBE8EcOk70DqHh3TgBAeRH8ABAxlHpQlDv6R/TIb1/XpLsazLT+sqXanGirdLMAFIHgR8Hu6B/Rw3uPTj+edNfDe4/q4b1HFedEMVAzKPWgYI/89vW8zzH8A1A7CH4UbHKOAf1SE5O6fcd+wh+ocgQ/CtZgc0/3PumujduHtaxnl268/4UytApAsQh+FGz9ZUuLWn/P4VOEP1CFCH4UbHOiTTdd3lLQkX/GnsOnmO8XqDKBTcQShM7OTh8cHKx0M1Cg/qGkbt+xf87avzQ1JtCWNW30+gFCUHMTsaB2JTriuucLl8w516/EuP9ANSH4UZLMhC+NBfyflBxPaXnPLko/QIUR/ChZoiOul79zrVZfdN6c6zLlI1B5BD8Cs+2WVXqt91ptXdc+Z/mH0g9QOQQ/Ajdzvt98KP0AlcFYPQhF9rj/q3t3K5lnXt/s0k/mdQDCxRE/Qtfd1VpQ6ef2Hfv5FwBQBqEf8ZvZa5L+LGlS0nvF9DVFfZg55WO+Xv+Z6wH4FwAQrnId8V/p7u2EfnQlOuLa03OVXu29VvGm2Jzrc/IXCA+lHpRdIaUfaerIf3Xvbso/QMDKEfwu6Wkz22dmG2Y+aWYbzGzQzAbHxsbK0BxU2sxeP/nG/jFNhX/mBPBt24d1R/9IOZsK1KXQx+oxsyXuftzMLpD0jKR/dvfnc63LWD3R1D+U1KadI0pNTE4vMynnuQCTdO+6dmr/QJaqG6vH3Y+n/56U9JikS8PeJmrLzH8BxJtieU8Au0TvH6BEoR7xm9lCSQvc/c/p+89Iutvdn8q1Pkf8yJit7382k3Tj5S1M+I5Iq7Yj/gsl/cbM9kv6naRd+UIfyNbd1TrrVb8ZLmnb3qMc+QNFCDX43f0Vd78kfVvp7t8Oc3uoH4mOuG68vKXg8KfrJ1A4unOiam1OtOnede1z9v6RpOMFlIUATGEGLtSM/qGkbts+nPPE76JzG+Uujacmph/f+fmV9P5BJFRbjR8ITL7yT2OD6f/enpgOfUl68+0Jbdw+rPZvPU39H5iB4EdNmVn+iTfFtPDss3Q6z/rjqQkmfQFmoNSDmre8Z1fefv/Z4k0xdXe1Uv5B3aHUg8hZUsCgbxLDPgAZBD9qXndXqxoXFNLxc6rr58N7j2p5zy5+ABBZBD9qXqIjrr61l6gp1ljwazI/ADfe/0J4DQOqFDV+1J3+oaT6BkYLGvJBousnal+xNX6CH3Vrtn7/uTQukPrWMvInag8nd4G0YoZ9kKSJ09LG7cP6SM8uun+irhH8qGuZfv8Lz557xq+M05r6ASD8Ua8IftS9REdcB+++RqsvOq+o123cPsyY/6hLBD8iY9stq7R1XXtRvX+S4ymu/EXdIfgRKYmOuIbv/ExRR/+piUlt5MIv1BGCH5G07ZZVuunylqJe8/Deo1r5zac4+kfNI/gRWZsTbXqt91qtuGBhwa95691Jhn1AzSP4EXnPfO0KbU2P+Clpzu6fmat+O+5myGfUJi7gAmboH0pq4/bhgtZlsndUg6q7gMvMrjGzUTM7ZGY9YW8PKFWiI15w/Z/J3lGLQg1+M2uQ9H1Jn5V0saT1ZnZxmNsEgrA50aabipjsfeP2YS3r2cWgb6gJYR/xXyrpkLu/4u7vSnpU0nUhbxMIxOZEm17tvbbgHwBJ2nP4FOGPqhd28MclvZ71+Fh62TQz22Bmg2Y2ODY2FnJzgOJlhn0o9MKvPYdPaXnPLq76RdUKO/hzHSi972yyu9/n7p3u3tnc3Bxyc4D5yVz4VUz5h6t+Ua3CDv5jkpZmPf6wpOMhbxMITfZk74VITUyqb2A05FYBxQk7+H8vaYWZLTezsyXdIOmJkLcJhCrREdeenqsKHvYhOZ6i9IOqEmrwu/t7km6VNCDpJUk73P1gmNsEymXbLasKDv9M6YerflENuIALCEj/UFKbdo4oNTE563om6d51zPSF4FTdBVxAVCQ64tqypk3xptisJ4Az/f4p/aBSOOIHQrK6d/ecE76bpn4I4k0xdXe18q8AzAtH/ECV6O5qLWjAN4munygvgh8ISbGTvdP1E+VCqQcIWf9QUn0Do3OWfWY6t3GBvrPmY5R/MCdKPUCVyfT737quXbHGhoJf9/bEaaZ8RCjOqnQDgKjIHLlnjv4zJ3bn8vDeo5LEmP8IDMEPlFGiIz79A5ApAR0fT835A7Bt71F1/s15lH0QCEo9QIVkSkCv9l6rBpv9FLBLnPhFYAh+oAqsv2zpnOskx1Na3bubcX9QMoIfqAKbE21zjvtjmgp/xv1BqQh+oEpsu2WVtq5r18Kzz+z5k+tEsGvqxG/H3U9z9I+i0I8fqELZJ36XNMUKugYg1rhAW+j3H0nF9uOnVw9QhbJ7/0iFjfuTmjit7p/un349kA+lHqAGFDLujyRNnHZt3D6s5T27qP8jL4IfqAHFjvuTqf9/+nu/CrFVqFUEP1AjMvP9NsUaC37Nyyff0spvPsXJX7wPwQ/UkERHXMN3fkY3Xd5S8GveeneSIZ/xPgQ/UIM2J9q0dV27Fp1b2NE/Qz4jW2jdOc3sLkm3SBpLL/q6uz8522vozgkUr38oqY3bhwtev8FM6y9byqBvdaTahmW+193b07dZQx/A/CQ64kWVfibdOfEbcZR6gDqQKf2c21j4V/rlk29pWc8uav8RFHap54uS/iRpUNLt7v5mjvU2SNogSS0tLR8/cuRIKO0BoqSYIZ8ladG5jbrz8yu58KtGFVvqKSn4zexZSR/K8dQ3JO2V9EdNdSn+N0mL3f2fZns/avxA8C7a9KQmC/iexxobtGVNG+Ffg8o6ZIO7X13IemZ2v6RflLItAPOz/rKl07N4zSY1ManbdzDkQxSEVuM3s8VZD6+XdCCsbQHIb3OiTSsuWFjQupPu9PmPgDBP7n7XzEbM7EVJV0q6LcRtAZjFM1+7Qhd+8OyC1k1NTGrj9mEme6ljDMsMRMiN97+gPYdPFbx+Zh6AeFNM3V2tlICqFMMyA8hr2y2rpu/3DyV1+479s574zTyTHE+p+2fU/+sF/fiBiEp0xHXPFy5RrPHMGb9ymZh0feu/DobcKpQDR/xAhGWO3vsGRgua5evNtyfCbhLKgCN+IOISHXHt6blKW9e1F3z0j9rGET8ASYUd/RczFwCqF8EPYFpmrt/+oaS6f7pfE6f/cuK3cYHprn9YOf145oTw9PqpHXTnBJDTbMHeP5TUpp0jSk1MnvG61Red977eQwhfWcfqCRrBD9SG1b27Zz0ZTPiXV7WNxw+gDh2fowdQMReJofwIfgBFW9IUq3QTUAKCH0DRurta6fpZwwh+AEVLdMS1ZU2bzjkrd4Ssvui8MrcIxSD4AcxLoiOu0c2fPSPkc53Y7R9KanXvbi3v2cWon1WAXj0AQpWv62escYG2rPkYff8DQK8eAFWlb2A0Z3//1MRpdf90P0f/FUDwAwjVbF0/J067+gZGy9gaSAQ/gJDN1fVzrmsCEDyCH0CourtaZbM8zzUB5UfwAwhVoiOuGy9vyflc4wJTd1drmVuEkoLfzNaa2UEzO21mnTOe22Rmh8xs1My6SmsmgFq2OdGmrevatejcvwzr3BRrVN/aS+jVUwGlDst8QNIaST/KXmhmF0u6QdJKSUskPWtmH3X3M0/tA4iEzJDPqLySgt/dX5IkszMqeNdJetTd35H0qpkdknSppBdK2R6A6GC8//CEVeOPS3o96/Gx9LIzmNkGMxs0s8GxsbGQmgOglmQu+kqOp+SSkuMpbdo5Qp//gMwZ/Gb2rJkdyHG7braX5ViW8xJhd7/P3TvdvbO5ubnQdgOoY7ku+kpNTNLnPyBzlnrc/ep5vO8xSUuzHn9Y0vF5vA+ACMrXtz85ntLq3t2Uf0oUVqnnCUk3mNk5ZrZc0gpJvwtpWwDqTL6+/SZR/glAqd05rzezY5JWSdplZgOS5O4HJe2Q9AdJT0n6Kj16ABQq13j/pjPrxZR/5qfUXj2PSXosz3PflvTtUt4fQDRlyjfZvXryzfHLkA/FK7UfPwCEYma//3wTvDPkQ/EYsgFATchV/ok1NjDkwzxwxA+gJuQq/9CrZ34IfgA1o5hhH7jyNz9KPQDqTq4rf2/bPqw7+kcq3bSqQPADqDu5rvx1Sdv2HqXfvwh+AHUoXxdPl+j3L4IfQB2arYsn/f4JfgB1aLbpHun3T/ADqEOZ6R5nhj/9/qcQ/ADq0uZEm+5d1654U0wmKd4U05Y1bXTpFP34AdQxpnvMjSN+AIgYgh8AIobgB4CIIfgBIGIIfgCIGIIfACKG4AeAiCl1svW1ZnbQzE6bWWfW8mVmljKz4fTth6U3FQAQhFIv4DogaY2kH+V47rC7t5f4/gCAgJUU/O7+kiSZ5RsOCQBQbcKs8S83syEz+7WZfSLfSma2wcwGzWxwbGwsxOYAAKQCjvjN7FlJH8rx1Dfc/fE8LzshqcXd/9fMPi6p38xWuvufZq7o7vdJuk+SOjs7vfCmAwDmY87gd/eri31Td39H0jvp+/vM7LCkj0oaLLqFAIBAhVLqMbNmM2tI3/+IpBWSXgljWwCA4pTanfN6MzsmaZWkXWY2kH7qk5JeNLP9kn4m6Svufqq0pgIAglBqr57HJD2WY/nPJf28lPcGAISDK3cBIGIIfgCIGKZeBIAK6h9Kqm9gVMfHU1rSFFN3V2vo00US/ABQIf1DSW3aOaLUxKQkKTme0qadI5IUavhT6gGACukbGJ0O/YzUxKT6BkZD3S7BDwAVcnw8VdTyoBD8AFAhS5piRS0PCsEPABXS3dWqWGPD+5bFGhvU3dUa6nY5uQsAFZI5gUuvHgCIkERHPPSgn4lSDwBEDMEPABFD8ANAxBD8ABAxBD8ARIy5V880t2Y2JunIPF9+vqQ/BticasA+1YZ63CepPverXvdpobs3F/qCqgr+UpjZoLt3VrodQWKfakM97pNUn/vFPk2h1AMAEUPwA0DE1FPw31fpBoSAfaoN9bhPUn3uF/ukOqrxAwAKU09H/ACAAhD8ABAxNR38ZrbWzA6a2Wkz68xavszMUmY2nL79sJLtLFa+/Uo/t8nMDpnZqJl1VaqNpTCzu8wsmfX5fK7SbZovM7sm/VkcMrOeSrcnCGb2mpmNpD+bwUq3Z77M7AEzO2lmB7KWnWdmz5jZy+m/iyrZxmLl2aeiv081HfySDkhaI+n5HM8ddvf29O0rZW5XqXLul5ldLOkGSSslXSPpP82s4cyX14R7sz6fJyvdmPlI/7f/vqTPSrpY0vr0Z1QPrkx/NrXc5/1BTX1PsvVIes7dV0h6Lv24ljyoM/dJKvL7VNPB7+4vuXu4sxJXwCz7dZ2kR939HXd/VdIhSZeWt3XIcqmkQ+7+iru/K+lRTX1GqALu/rykUzMWXyfpofT9hyQlytqoEuXZp6LVdPDPYbmZDZnZr83sE5VuTEDikl7PenwsvawW3WpmL6b/6VpT/9zOUk+fRzaX9LSZ7TOzDZVuTMAudPcTkpT+e0GF2xOUor5PVR/8ZvasmR3IcZvtyOqEpBZ375D0NUk/MbO/Lk+LCzPP/bIcy6qyP+4c+/cDSRdJatfUZ3VPRRs7fzXzeRRptbv/vaZKWF81s09WukGYVdHfp6qfetHdr57Ha96R9E76/j4zOyzpo5Kq5kTVfPZLU0eUS7Mef1jS8WBaFKxC98/M7pf0i5CbE5aa+TyK4e7H039Pmtljmipp5TqPVoveMLPF7n7CzBZLOlnpBpXK3d/I3C/0+1T1R/zzYWbNmZOeZvYRSSskvVLZVgXiCUk3mNk5ZrZcU/v1uwq3qWjpL1zG9Zo6mV2Lfi9phZktN7OzNXXi/YkKt6kkZrbQzD6YuS/pM6rdzyeXJyTdnL5/s6THK9iWQMzn+1T1R/yzMbPrJf27pGZJu8xs2N27JH1S0t1m9p6kSUlfcfeST4iUS779cveDZrZD0h8kvSfpq+4+Wcm2ztN3zaxdU2WR1yR9ubLNmR93f8/MbpU0IKlB0gPufrDCzSrVhZIeMzNpKh9+4u5PVbZJ82Nmj0i6QtL5ZnZM0p2SeiXtMLMvSToqaW3lWli8PPt0RbHfJ4ZsAICIqctSDwAgP4IfACKG4AeAiCH4ASBiCH4AiBiCHwAihuAHgIj5f8PKw2fPl/gPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X2[:,0],X2[:,1])     # 第二个维度上的主成分分布在一条直线上，这是一种极端情况，最终求出的主成分也会坐落于这条直线上\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.62997951,  0.77661175])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 求第二个维度上的主成分\n",
    "w2 = first_component(X2, initial_w, eta)\n",
    "w2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.14402316262219e-06"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 验证w和w2是否垂直，可使用两者相互点乘（点乘即两者的模相乘，再乘以两者夹角的余弦，且知cos90度=0）\n",
    "w.dot(w2)  #结果非常接近于零，排除掉计算机计算浮点数的误差，可以认定w和w2是互相垂直的"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 整合PCA算法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def first_n_components(n, X, eta=0.01, n_iters=1e4, epsilon=1e-8):\n",
    "    \n",
    "    X_pca = X.copy()\n",
    "    X_pca = demean(X_pca)\n",
    "    \n",
    "    res = []\n",
    "    \n",
    "    # 求前n个主成分\n",
    "    for i in range(n):\n",
    "        # 随机生成一个初始化的w\n",
    "        initial_w = np.random.random(X_pca.shape[1])\n",
    "        w = first_component(X_pca, initial_w, eta)\n",
    "        res.append(w)\n",
    "        \n",
    "        # 每计算出一个维度上的主成分，就将其减去再继续求下一维度的主成分\n",
    "        X_pca = X_pca - X_pca.dot(w).reshape(-1,1) * w\n",
    "        \n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([0.77660931, 0.62998252]), array([-0.62997973,  0.77661157])]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "first_n_components(2, X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
